Steady state particle distribution of a dilute sedimenting suspension 
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Sedimentation of a non-Brownian suspension of hard particles is studied. It is shown that in 
the low concentration limit a two-particle distribution function ensuring finite particle correlation 
length can be found and explicitly calculated. The sedimentation coefficient is computed. Results 
are compared with experiment. 



I. INTRODUCTION 

In recent years the familiar problem of sedimentation 
of a non-Brownian suspension of particles has gained in- 
terest as new insight into the phenomenon is gathered. 
Despite vivid research in the field, some questions still 
remain intriguing puzzles. Even when considering obser- 
vation time scales and dimension regimes, which sanction 
the neglect of Brownian motion, the dynamics show to 
be many-body and highly chaotic leading to astonishing 
consequences. 

One of the main, still open problems, is the determi- 
nation of the particle probability distribution function 
for a steadily sedimenting suspension. It was first chal- 
lenged by Batchelor I*!. In order to derive an equation 
for the pair-distribution function he considered a polidis- 
perse suspension. Taking then the limit of identical par- 
ticles, a particle distribution function for mono-disperse 
suspensions could be derived. But these results where 
shown to be ambiguous (depending on the way the limit 
is taken) and therefore doubtful (Q, 0, |1]). The prob- 
lem of particle distribution is often skimmed over by the 
assumption, that the particles are randomly distributed 
in the fluid. Indeed this will be the case, if Brownian 
motions play a significant role in the dynamics of the sys- 
tem, creating a randomising mechanism. But whenever 
this type of stochastic motion is absent, there seems to 
be no justification for such assumptions. Further, it has 
been '4^ pointed out, that a equilibrium (equal probabil- 
ity) distribution of particles results in characteristics of a 
non-Brownian suspension, which have not been detected 
in experiments (51 . Moreover it can actually be explicitly 
shown that such a distribution cannot be a solution of the 
Liouville equation incorporating full multi-body interac- 
tions, therefore falling completely out of consideration. 

In this work we propose a way to tackle this issue. The 
derivation of the statistical properties of the system will 
be based on the conjecture that correlations in steady 
state must be of finite length. This assumption, orig- 
inating in the basic ideas of statistical mechanics, has 
profound consequences. It, as will be shown, leads to 
an effective equation for the pair distribution function. 
Truncation of long range hydrodynamic interactions will 
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be achieved through a screening mechanism similar in 
spirit to that introduced by Koch&Shaqfeh 0. 



II. BBGKY HIERARCHY 

The system under consideration is a suspension of N 
spheres of equal radius a immersed in an incompressible 
fluid of shear viscosity 77. In the regime studied the par- 
ticle Reynolds number is assumed to be small enough to 
treat inertial effects as virtually absent. Moreover the 
Peclet number describing the relative influence of hydro- 
dynamic effects to Brownian motion, is presumed to be 
large. On those assumptions and on the laboratory time 
scales imposed the fluid motion is governed by the sta- 
tionary Stokes equation. The dynamics of the particles 
are therefore given by the equations of motion 

-^ = U^ = Y.^^.AX)F, (1) 

i=i 

where fj. is the translational part of the mobility ma- 
trix, principally a matrix function of the configuration 
X = {Ri, . . . , Rn) of all the particles (if not indicated 
otherwise, dimensionless - normalised to 2a - distances 
will be used). It describes the hydrodynamic interac- 
tions between particles and connects the forces acting on 
particles to the velocities they acquire in a given config- 
uration. The forces Fj will be assumed to be constant 
and equal F for all particles. 

The particle distribution function satisfies the Liou- 
ville equation, which for the system considered, has the 
following structure 

+ ^ V. . {X)Fp{X; t)] = (2) 

The mobility matrix can be found following a formal pro- 
cedure presented for example in rcf. 7J. Summarising, 
in order to calculate the hydrodynamic interaction, the 
boundary conditions on the particles surfaces are sub- 
stituted by induced force densities introduced for each 
particle. These force densities, induced by fluid flow, 
produce subsequent fluid velocity fields which are prop- 
agated and reflect off other particles and thereby cause 
interactions. In consequence, the mobility matrix can 
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be expressed as a scattering sequence equivalent to a su- 
perposition of all these reflections. This sequence con- 
tains one-particle scattering operators, which construct 
induced force densities on particles given velocity fields 
around them, and Green operators which propagate the 
influence of a force density concentrated on a given par- 
ticle, resulting in fluid velocity fields around other parti- 
cles. In an unbounded, infinite system, the Green opera- 
tor is the Oseen tensor. 

The full TV-body mobility matrix can therefore be ex- 
panded in terms of the number of particles which enter 
each term of the scattering sequence. This so called clus- 
ter expansion has the structure 



N 



i=2 
N 



/Xi2(X) = /xp2'(12) + ^Mf2(12i) 



(3) 



j=3 



where (1 . . . , s) denotes all the terms of the scattering 
sequence of /x which involve only, but all the particles 

Next we introduce reduced distribution functions 



n(l,2,...,s) = 



{N~sy. 



d{s + !)■■■ dNp{X) (4) 



These functions represent the average densities of pairs, 
triplets, etc. of particles in given configurations. 

An integration of the Liouville equation in all the vari- 
ables except the first s, together with the cluster expan- 
sion of the hydrodynamic interactions, leads to an infi- 
nite hierarchy of equations governing the time evolution 
of the reduced distribution functions. This is the ana- 
logue of the BBGKY (Bogolubov-Born-Green-Kirkwood- 
Yvon) hierarchy introduced in the kinetic theory of gases 
It should be considered in the thermodynamic limit 
when the number of particles and the volume of the sys- 
tem go to infinity, while the density is kept constant. 
The equations for the two and three particle distribution 
functions have the structure 



dn{l2) 
dt 



V,-/x,^.(12)Fn(12) 



[ d-iY^V,- [/xF],(12; 3)n(123) + . . (5) 



i=l,2 



where for example 

[mF]i(12;3) = [m1i^(13) + m13^(13) + 5]m1?(123)]F (6) 



and 



an(123) 



= - ^ V,-M,^.(123)Fn(123) 



(7) 



where in both equations, terms depending on higher- 
order distribution functions have been omitted. 

The reduced distribution functions factorize for groups 
of particles that are far away from each other. Therefore 
these functions can be expanded in terms of correlation 
functions according to 



n(12) = h{l)h{2) + h{\2) 
n(123) = h{l)h{2)h{2,) + h{l2)h{i) 

+h{li)h{2) + h{2'd)h{l) + /i(123) 



(8) 



where . . . s) is a s-particle correlation function which 
vanishes whenever any subset of particles is dragged away 
from the rest. If the correlations in a system are to be 
of finite length, which will be assumed, the correlation 
functions must decay faster then the inverse of the inter- 
particle distance cubed. 

Using the expansion (|SJ|, the BBGKY hierarchy for the 
reduced distribution functions can be transcribed into a 
hierarchy of equations for the time evolution of the suc- 
cessive correlation functions. An analysis of these equa- 
tions shows that if system configuration is such that par- 
ticles form two uncorrelated clusters, some expressions on 
the r.h.s, formulated in terms of the scattering sequence, 
contain single Green operators (so called connection or 
articulation lines 0j[3) joining these two groups of par- 
ticles. Such terms result in long range contributions be- 
cause a solitary Green operator includes slowly decaying 
parts proportional to where r is the relative dis- 

tance between the groups and 7 = 1,2, 3. The evolution 
of a correlation function is therefore given by an equation 
which contains long range terms and is therefore non- 
local in space. This structure of the hierarchy equations 
is inconsistent with the finite correlation length assumed 
for the correlation functions. As will be shown in the 
next section, considering the low concentration limit a 
particle distribution can be found, which cancels out all 
long range terms in the hierarchy, thereby saving finite 
correlation length. 



III. LOW CONCENTRATION LIMIT PAIR 
DISTRIBUTION FUNCTION 

In the low concentration limit the dominating part of a 
s-particle reduced distribution function is proportional to 

. Consequently, the first terms of the hierarchy equa- 
tions have a well established order in density if the con- 
centration of particles is low. 

Consider the equation giving the time evolution of the 
two-particle correlation function. It is derived through a 
substitution of © into ||SJ). In the lowest order of density 
it vanishes due to the relation Q 



Y V,-/x,^(12)=0. 



(9) 



which is a consequence of the isotropic nature of the 
hydrodynamic interactions. The two-particle equation 
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yields therefore in the dilute limit no condition for the 
two-particle correlation function. This peculiarity was 
also encountered by Bathelor . His solution was based 
on the analysis of a polidisperse suspension where the 
two-particle equation yields a condition, but was shown 
to be ambiguous (0 We turn our attention to 

the next hierarchy equation, i.e. the third equation which 
emerges when ^ is addressed to {TJ, truncated to the 
terms proportional to the density cubed. 

As pointed out, the r.h.s of this equation contains long 
range terms. If all three particles are far away from each 
other, both sides of the equation contain only shot range 
terms. But consider the case when one particle (e.g. the 
particle with index 1) is far away from a close pair of 
particles (consequently particles with indices 2 and 3). 
All the terms, which in such a configuration lead to long 
range contributions can be identified and singled out re- 
sulting in the expression 



nV2 ■ [(At2i(2|l) 
nVa- [(M3i(3|1) 



/X2i(23|l)-f /X2i(2|3|l)M23)] + 
/X3i(32|l)+/X3i(3|2|l))n(23)] (10) 



where each vertical line stands for a single Green op- 
erator connecting uncorrelated particles. For example 
/X2i(23|l) stands for the sum of all terms of the scatter- 

ing sequence of (123) that contain only a single Green 
operator connecting particle 1 and the group (2,3). 

In what follows, we will show that there exists a two- 
particle probability distribution function which cancels 
out these long range terms and therefore leads to finite 
correlation length 

In equation H10|l the long range contributions appear 
in the Green operators binding particle 1 with the group 
consisting of particles 2 and 3. Their explicit form is 
retrieved through an analysis of the multipole matrix el- 
ements of the Green operators, which can be found e. g. 

in g]. 

When taking into account the symmetry relations 
upon exchange of particles, and a change of variables 
where r becomes the vector joining the distant particle 
and the center of the group of close particles, while R 
the vector denoting the relative position of the two par- 
ticles within the group, we reach the conclusion that most 
of the long range contributions cancel out automatically 
- only the multipoles proportional to the inverse of the 
inter-particle distance survive. 

A further expansion in terms of i? = |i2| {R ^ r) 
shows that only terms proportional to 1 /r^ remain. The 
equation which emerges has an identical structure to the 
one derived by Batchelor&Green and describing the 
relative motion of a pair of particles in a shear flow. The 
resulting particle distribution function n{R) = •n?g{R) is 
isotropic. It is given by the solution of the differential 
equation 



{I -A) 



1 dg i{A - B) 
'g'dR 



R 



dA 
dR 



= 



(11) 



A and B are functions of the scalar distance R given by 
the relations A = 2{a\i+a{fj/R and B = 2{(3\i+ I3lf} / R 



Figure 1: The two-particle reduced distribution function g{R) 
satisfying the finite correlation length criterion 
9{R) 



R 



where a*^ and /3*^ are functions appearing in the explicit 
form of the generalised two-particle mobility matrix con- 
necting the translational velocity of the particles to the 
stresslet Since Vmvji^oo 9{R) — 1 



9{R) = 



1-A 



cxp 



3(B ~ A) 
R R{1 - A) 



dR 



(12) 



Note that this solution is independent of the direction dis- 
tinguished by the external force field, the effects of which 
are in the low concentration limit completely dispersed 
by the isotropic hydrodynamic interactions. Further, to 
verify the consistency of the reasoning imposed, it can be 
shown that the two-particle distribution function derived 
(|12|l guarantees convergence of all terms up to order 
in equation (jSJ governing the time evolution of the two- 
particle reduced distribution function. 



IV. RESULTS 

A. Stationary particle distribution function 

The functions A{R) and B{R) can be expressed as 
series in inverse powers of inter-particle distance |llj |. 
In numerical approximations of the integral (|12|l these 
sequences are truncated at 1000 and 800 terms respec- 
tively. Lubrication corrections and far limit asymptotic 
with g{R) — 1 « 0.1953/i?^ are taken into account. The 
computed function g{R) is plotted at fig.d Further, the 
structure factor at fc = is calculated: S{0) = 1 — 1.640. 



B. Sedimentation coefficient 

In the low concentration limit the average sedimenta- 
tion velocity [/ of a particle in the suspension measured 
relative to the Stokes velocity Uq may be expanded in a 
series in the powers of volume fraction (f) = Ana^n/S. 



^ = l + 0/f + O((/.2) 



(13) 
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Figure 2: Dimensionless sedimentation velocity as a function 
of volume fraction - experimental results after Hanratty et 
al.|l4ll. Points represent experimental data and the solid curve 
is a fit. The solid line represents the sedimentation coefficient 
K = —3.87, whereas the dotted line is the coefficient for a 
suspension in equilibrium with K = —6.55. 
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The linear coeffieient K ean be expressed in terms of a 
microscopic expression 13] involving two-body mobility 
matrices and the pair distribution function. 



K = Ka + / (.g(12) - l)Tr 



1=1 



dRi2 
(14) 

where Kq = -6.546 (Jl^, fl2() is the sedimentation coef- 
ficient for a equilibrium distribution of hard spheres and 



fiQ is the mobility of a single sphere. Dimensionless dis- 
tance Ri2 normalised to 2a is adopted. The integral can 
be interpreted as a correction emerging due to the change 
of the distribution away from equilibrium. A numerical 
calculation yields the result X = — 3.87. It is compared 
with experimental data of Hanratty et al. on fig. El A 
fit to the experimental data of Ham et al. leads to 
the sedimentation coefficient equal to —3.9, in excellent 
agreement with the calculated value. 



V. CONCLUDING REMARKS 



We have shown that a mono-disperse suspension 
of hard, non-Brownian particles can develop a micro- 
structure which insures finite correlation length. A 
scheme based on the formal schemes of non-equilibrium 
statistical mechanics and a detailed analysis of hydro- 
dynamic interaction was used in the low-concentration 
regime to derive an equation for the two-particle distri- 
bution function. The solution shows to be isotropic due 
to the domination of the isotropic hydrodynamic inter- 
actions. The mechanism which insures the cutoff of long 
range correlations can be interpreted as hydrodynamic 
screening in the sense that the micro-structure of the 
suspension arranges itself in a way which truncates long 
range parts of the interactions. A close pair of particles 
feels the influence of a third distant particle through an 
effective shearing flow. Screening leads to a conflguration 
of the close pair, which hinders the effect of the distant 
particle. The idea, that a suspension might develop a mi- 
crostucture which changes the characteristics of it, was 
introduced by Koch&Shaqfeh 651 upon the analysis of the 
problem of diverging velocity fluctuation. 

Further, the sedimentation coefficient was calculated. 
The value —3.87 was found. It agrees very well with the 
experimental results cited. 
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